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Abstract 

Dynamics of individual DNA undergoing constant field gel electrophoresis (CFGE) is studied by 
a Brownian dynamics (BD) simulation method which we have developed. The method simulates 
electrophoresis of DNA in a 3 dimensional (3D) space by a chain of electrolyte beads of hard spheres. 
Under the constraint that the separation of each pair of bonded beads is restricted to be less than 
a certain fixed value, as well as with the excluded volume effect, the Langevin equation of motion 
for the beads is solved by means of the Lagrangian multiplier method. The resultant mobilities, 
[/,, as a function of the electric field coincide satisfactorily with the corresponding experimental 
results, once the time, the length and the field of the simulation are properly scaled. In relatively 
strong fields quasi-periodic behavior is found in the chain dynamics, and is examined through the 
time evolution of the radius of the longer principal axis, Ri(t). It is found that the mean width of 
a peak in Ri(t), or a period of one elongation-contraction process of the chain, is proportional to 
the number of beads in the chain, M, while the mean period between two such adjacent peaks is 
proportional to M° for large M. These results, combined with the observation that the chain moves 
to the field direction by the distance proportional to M in each elongation-contraction motion, yield 
\x oc M°. This explains why CFGE cannot separate DNA according to their size L (ex M) for large 
L. 
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I. INTRODUCTION 



Gel electrophoresis is a most widely used technique to separate DNA according to their 
length L under an electric field.S The constant field method is quite efficient for DNA of 
small sizes. It is well known, however, that the separation does not work well for DNA of 
sizes longer than a certain value which depends on the field and the "pore" diameter of the 
gel.i The non-constant field methods such as the pulsed-field gel electrophoresisi have been 
empirically invented and they have remarkably improved the shortcoming. Nevertheless, 
the understandings on these gel electrophoresis methods, even on the constant field method, 
still remain unsatisfactory from the view point of statistical physics of polymers. 

By constant field gel electrophoresis (CFGE), DNA molecules with the radius of inertia, 
Ri, smaller than the pore diameter of gel, &gel, £1X6 sieved by gel.BS In this sieving regime DNA 
with the smaller Ri migrate the faster through the gel. In the opposite case with Ri > a ge i, 
DNA is considered to move through a "primitive chain"! which is determined by stochastic 
movement of its front end in pores of the gel under the field. Si'i This picture, denoted as 
the biased reptation model, successfully describes various features of CFGE, such as the 
L-dependence of the migration mobility, fi, observed experiment ally@ in sufficiently weak 
fields. More recently, the model of the biased reptation with fluctuations has been proposed0 
which well describes the observed field dependence of \x in moderate fields.0 When the field 
is further increased, the field-dependence of \i begins to be saturated. In this regime, inter- 
esting behavior of an individual DNA has been observed: "hernias" -like configurations by 
the computer simulations0'0 and "quasi-periodic" behavior, i.e., "elongation-contraction" 
motion in a cyclic way by the real-time fluorescent microscopy experiments.01l! Indeed, 
the latter experimental method combined with optical or magnetic tweezers for a molecule 
reveals that the elasticity of DNA is originated from the entropy associated with its static 
conf or mat ions However, dynamic properties of DNA under CFGE are not described 

in a unified way by this picture alone. 

In order to get further insights into the problem mentioned above, we have developed a 
Brownian dynamics (BD) method which simulates electrophoresis of DNA in a 3 dimensional 
(3D) space. The method has been extensively applied to the CFGE, in particular, the 
elongation-contraction motion of the chain under CFGE, the results of which we will discuss 
in the present series of papers. In the present paper I our numerical method is explained and 
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global aspects, such as the mean mobility \i and time evolutions of the velocity of the center 
of mass, vc(t), and of the radius of the longer principal axis, Ri(t), are investigated. In the 
accompanying paper we analyze dynamics of 'defects' introduced by de GennesEI in a 
chain under CFGE. Finally in paper hi! we will examine conformational changes of the 
coarse-grained chain under CFGE in details, and will argue that the field-parallel component 
of its elongation-contraction motion is essentially understood as a deterministic dynamics of 
an elastic string in a ID space with an obstacle, thereby an origin of the elasticity will be 
shown to be the conformational entropy of the chain. 

In the Brownian dynamics (BD) simulations extensively performed in the present series 
of studies, DNA is modeled as a chain of spherical electrolyte beads with a constant ra- 
dius 1/2 in a 3D space filled with a solvent. Each bead interacts with other beads of the 
chain and spatially-fixed beads constituting gel by hard-core repulsive interactions. It also 
interacts with its two bonded beads by a non-linear elastic force which keeps the distance 
between the two to be less than a threshold value which we put v2- To solve the Langevin 
equation of motion under these constraints as well as under the field and random forces, 
we employ a BD algorithm with the Lagrangian multiplier method. It is similar to the one 
by Deutsch and Madden,0El but is essentially different from their method in the following 
respects. In their method, the time pitch to solve the equation is adjusted so as to keep 
a fixed distance between bonded beads under the tensile forces which are evaluated by the 
Lagrangian multiplier method. In our method, on the other hand, the time pitch is fixed 
at a certain value, which we regard as the mean period of random forces. Within each time 
step we evaluate paths of all beads affected by the constraints, i.e., by hard-core interac- 
tions and extremely nonlinear elastic interactions between a pair of beads. These effects are 
reread as the supposed constraining forces with the appropriate Lagrangian multipliers in 
an iterative way. The resultant constraining forces thus specified yield a new set of bead 
positions without any violation to the constraints. To our knowledge the present BD anal- 
ysis is the first extensive simulation of DNA gel electrophoresis in a 3D space based on the 
Langevin equation involving purely mechanical, microscopic forces which are supposed to 
act on monomers of DNA. 

The purposes of the present paper I are to explain our BD method and to demonstrate 
its results on some global properties of DNA under CFGE. The resultant mobility fi as a 
function of the field turns out to agree well with those observed experimentally once a set of 
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three parameters, /*, r* and E*, are properly specified. They relate respectively the scales 
of length, time and field strength of the simulation and the experiments. 

In the BD simulation under relatively large fields quasi-periodic time evolution is observed 
in Ri(t) (see Fig. |5| below), as well as in vc(t), as has been observed experimentally.B'El Each 
local maximum of Ri(t) at t™ ax with n = 1,2, ... corresponds to a conformation that a chain 
just gets rid of trapping due to gel. The local minimum at t™ in:1 preceding this maximum is 
then considered to be the instance that, being trapped by gel, the chain starts to elongate 
from a certain coiled conformation, while at time £™ in . r of the local minimum subsequent 
to i™ ax the chain retrieves another coiled state. We establish a method to properly classify 
a whole time sequence of Ri(t) into peak parts, from £™ irrl to t™ in:r and the rest. We call 
the former parts deterministic ranges and the latter stochastic ones. We then investigate 
statistics of periods of the deterministic ranges, D n = t 7 ^^— t™ in:1 , and those of the stochastic 
ranges, (1/A) n = £^f^ a — £™ in:r . It is found that {D n } obey a Gaussian distribution with the 
mean D = D n oc M and the width A oc D, while {(1/A) n } do a Poisson distribution with 



the mean 1/A = (1/ A) n oc M° for large M, where M is the number of beads of the chain and 
the overline stands for the averages over n, i.e., peak structures. These results, combined 
with the observation that in each elongation-contraction motion the chain moves to the field 
direction by the distance proportional to M, yield fi oc M° and explain why CFGE cannot 
separate DNA according to their sizes L (oc M) for large L. 

The organization of the present paper is as follows. In the next section we explain our BD 
method. The mobility \i of the chain obtained by the BD method is presented and compared 



with the experimental results in Sec. |T|. The method of analysis of peak structures in Ri(t) 
is explained, and statistical nature of the deterministic and stochastic ranges is discussed in 
Sec. [TV]. The last section is devoted to further discussions of our results. 



II. THE BROWNIAN DYNAMICS METHOD 



A chain of spherical beads in a continuous 3D space, shown schematically in Fig. [I], 
is investigated by a new type of Brownian dynamics (BD) calculation. The hard core 
interaction is assumed in order for centers of any pair of beads not to come closer than 1 
(the excluded volume effect), where the radius of beads is set 1/2. Furthermore distances 
between neighboring beads, I, are restricted to I < We suppose that this is realized by 
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the extremely nonlinear elastic interaction between the neighboring beads: it is infinite for 
/ > \/2 and zero otherwise. These interactions assure the self-avoiding walk of the chain. 

The BD method we adopt in the present work is to solve the following Langevin equation 
of motion for the chain, 

(x t = qE h + J2J2 5 S + E I*? " ^-i] + C 1 ) 

a j a 

where £Cj is the position of the center of the z-th bead, and (, q, and E h are the viscosity, 
the charge of a bead and the (bare) external field, respectively. The second and the third 
terms in r.h.s. are the constraining forces required for the beads to satisfy the conditions, 

ll^ijll = ^* ^ (^) 

pill = ||aj i+ i - Xi\\ < V2. (3) 

The vectors Sfj and are respectively written as &^[xi — Xj] and <fif[xi — x i+1 ] by intro- 
ducing the Lagrangian multipliers 9fj and (f)f. The last term in eqn. [I] stands for the random 
force from a solvent. Its distribution is characterized by the following averages 

(fi) = (4) 
(fMfjit')) = 2(k B TAt6 ij 5{t-t'). (5) 

Here, fee, T, and At are the Boltzmann constant, the absolute temperature, and the average 
period of the random force, respectively. 

To solve the one step evolution of eqn. [I], thereby the terms {S°!j , F?} with the Lagrangian 
multipliers are determined, we examine the full Langevin equation with the inertia term 
and the forces arising from the hard core interaction and the nonlinear elastic interaction 
mentioned above. We solve the equation in a time interval of At > t > 0, and take the 
overdamped limit, m/( 0, where m is the supposed mass of the bead. The constraints 
of eqns. @ and § are supposed to be satisfied by positions of the beads at t — 0, denoted 
by {x®}. At this instance fi and qE^At are added as impulsive forces with the integrated 
strength fi where fi = fi + qE h At. When the constraints are discarded, the solution in the 
overdamped limit is described by the relation between the displacement Axi in this time 
interval and the force as 

Axi = jf u (6) 
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which is just what is obtained by simply integrating eqn. [I] without the second and the third 
terms. 

In the case that positions of a pair of beads x® 1 = x® + Axj and x ^ = x® + Axj 
violate the constraint of eqn. |2|, we solve the full Langevin equation for this pair with 
the hard core scattering, which can be either elastic or inelastic. We then obtain, again 
in the overdamped limit, positions xj and ccj at t = At as shown in Fig. in case of 
the inelastic scattering. The result is reread as the consequence of the constraining force, 
S}j = 0\a[x\ — x^], which is determined by the conditions x\ — = (/, ; + Sjj)/( and 
x j — x® = (fj — Sjj)/( (note that Sji = —Sh). This constraining force Sh just corresponds 
to the one in the second term of eqn. [I]. In the case x® 1 and ar^K violate eqn. ^|, the nonlinear 
elastic force is supposed to reduce the distance ||Zj|| to \/2. This effect is represented by the 
constraining force (= —Fl +1 ), which is added to the third term of eqn. [l]. Its Lagrangian 
multiplier <pj is determined in such a way that \\x\ — x\ +l \\ = \[2 with x\ = x® + F^At and 
xj +1 = x Q i+l + Ff +l At (see Fig. ||b). If x® 1 involves violation to more than one constraints of 
eqns. ^| and/or |3], the above-mentioned procedures are applied to evaluate the constraining 
forces for each pair of beads with the common x® 1 . If, on the other hand, x® 1 does not 
involve any violation to eqns. ^| and |3|, we put xj = x® 1 . 

In the configuration {xj}, it is generally expected that new pairs of beads may violate 
the constraints since we have tried to remedy the violations in {a;? 1 } in a pairwise way. 
Then we repeat the above-mentioned procedure, in which fi is replaced by f i + C\, where 
C\ = J2j Sjj + Fi ~ -^i-i 1S ^e constraining force determined by the first procedure. If a 
configuration without any violation to the constraints is obtained by the a con times repetition 
of this procedure, we regard {a?" con } as the solution of eqn. [I] evolved from {x®} by one unit of 
time At. In this manner the multi-scattering processes in the interval of At, including those 
due to the nonlinear elastic interaction, can be taken into account, and the corresponding 
forces {S°j, F^} in eqn. [I] are determined. 

Although our BD method is based on the overdamped Langevin equation with the La- 
grangian multipliers, it distinctly differs from the one due to Deutsch and Madden .IllJll In 
the latter the equidistant condition on all the pairs of beads is imposed, and the constrain- 
ing tensile forces associated with it are evaluated at the instance when random forces are 
applied. Then the time interval At is adjusted so that the configuration at t + At in fact 
obeys the equidistant condition within accuracy of 0.1%. In contrast, our BD method takes 
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into account a similar, but looser, constraint on the bond length as well as the excluded vol- 
ume effect of the beads by explicitly examining the supposed multiple scattering within the 
interval At which is fixed. Therefore there is no finer time unit other than At in solutions 
of our BD method. It is the mean period of random forces {/«}. 

In the present work we simulate dynamics of a chain in a continuous 3D space confined 
into a cubic box with a volume 360 3 . The periodic boundary condition is imposed to each 
direction of the box. The gel is represented by a network of immobile bars arranged so 
that they form a simple cubic lattice. Each bar consists of tightly connected beads with a 
radius 1/2, and its interaction with a bead of the mobile chain is computed according to 
the constraint of eqn. |2|. We employed three lattices of gel with different lattice distances 
a ge i = 10, 17, and 20, which are considered to correspond to the 'pore' diameter of agarose 
gel.i The lengths of chains, M, used in the BD calculation are M = 30, 40, 60, 80, 160, 
240, and 320. We have performed typically 64 runs whose initial configurations are given 
independently. The direction of the field is fixed at (1,1,1). The time is measured in unit 
of At, i.e., 1 BD step. 

Here we want to emphasize that the present BD method correctly reproduces the expected 
static and dynamic properties of a real (not phantom) polymer in a 3D space under the 
vanishing field: Ri oc M VF with vp ~ 0.59 ± 0.01 and the Einstein relation Dq oc /i/M are 
confirmed (see Fig. ^ in Sec. pTT| ) . Here, vp, Rj, and D G are the Flory exponent, the radius 
of inertia, and the diffusion constant of the center of mass, respectively. 

The number of the multiple scattering, a con above mentioned, is expected to be at most 
of the order of M 2 . It has turned out that a con becomes relatively larger when the chain is 
in an extended configuration under gel electrophoresis. We quote here typical figures for a 
chain with M = 80: a con ~ 4 x 10 1 , 5 x 10 2 for E = 0, 0.032, respectively. 

III. MOBILITY 

The most basic quantity that has been commonly measured in CFGE of DNA is the long 
time average of velocity divided by the electric field, or the mobility fi. In the present work, 
we have implemented long time BD simulations with a highly efficient algorithm using a 
massive parallel computer and have eventually obtained /i, hereafter denoted by /xbd ; with 
enough statistical precision for all the set of parameters investigated. Here we only mention 
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typical figures of computation to get /zbd of a M = 240 chain for each run: more than 
150 hours are consumed on a single processor of SGI2800 system to pass 10 6 BD steps. In 
the present work the following values of the parameters are chosen: At = 1, ( = 10 and 
k^T = 2. The results are shown in Fig. [5] with the experimental data due to Heller et a/.il 
Here and hereafter the value of field E stands for qE^ of eqn. [I]. 

The simulated hbd at E = in Fig. [3] is evaluated through the Einstein relation 



= Dg = } im h^Ro® ~ ^g(o)) 2 )u =0 , (7) 

M t^oo <ot 

where Rc{t) is the position of the center of mass of the chain. The diffusion constant 
Dq in the r.h.s. is evaluated by the BD simulation in a vanishing field in the time range 
t = (5 ~ 10) x 10 6 . The resultant /i BD turns out to be compatible with the value of 
liniE^o /ibd as seen in the figure. This confirms that the Einstein relation in fact holds in 
the present BD simulation for the chain. 

In order to compare the BD results with the experimentally observed data, we have to 
specify proper conversion of values of the parameters in the simulation and the experiment. 
Let us first consider the length scale of a chain. In the accompanying paper II we examine 
in detail the embodiment of 'defects' due to de Gennes0 in our BD chain and obtain the 
mean distance (n) of adjacent 'defects' along the chain under E = as (n) ~ 5.7. This 
value, in unit of the mean distance between neighboring beads, (I), is almost independent 
of the set of parameters a ge i and M we have examined. Regarding (n) as an estimate of the 
'persistent length' in the BD simulation, Vbd, we put Vbd = c(n), where c is an adjustable 
parameter whose value is nearly equal to unity (note that Vz > (/) > 1). The conversion 
factor of the length-scale I* is then determined by 

Pexp = 1*Pbd = cl*{n), (8) 



where V exp is the persistence length experimentally measured. With V exp ~ 60nm at T exp ~ 
300K taken from Ref. ^4], we obtain cl* ~ lOnm. 

Another length scale of importance is the mean pore size of gel, which is rather difficult 
to be accurately specified even in the experiments^ Here we simply suppose that DNA in 
1% agarose corresponds to a BD chain in the gel with a ge \ = 17. With Vbd above fixed, this 
value of a ge i is about three times larger than that of Vbd- 

Next let us relate E in the BD simulation to the experimental field, E exp , in V/cm. For 
this purpose we impose that the ratios As/ZcbT in the BD calculation and the experiment 
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be identical, where Ae is the energy that a part of chain of a length equal to the persistence 
length V gains when it moves a distance of V under the electric field: 



^B^exp &B^BD 

or 



(9) 



£ exp = E*E; E* = fcBTc *pfBD 5 1q2 -2 v/cm _ (1Q) 

where a is the density of charge of DNA, and /cr^bd = 2 has been used. Quoting further 
the experimental value of a ~ O.6e~/A,0 as well as cl* ~ lOnm at T exp ~ 300K in eqn. |8|, 
we obtain the last expression for E* in the above equation. As for the ordinates of Fig. |3| 
we may relate //bd to /i eX p[cm 2 /V • s] by 

/^ex P = /i> BD ; li* = ^ 2 x 10- 9 4 [cm 2 /V • s], (11) 

where r*[s] is the real time which corresponds to one discretized time step At in the BD 
calculation with £ = 10 (or At/( = 0.1). 

The two sets of data of /i in Fig. |3| are drawn as follows. The data /i exp are plotted directly 
using the units indicated on the right ordinate and the upper abscissa, whereas the units 
for /iBD are indicated on the left ordinate and the lower abscissa with c 2 = 1.60 and 1.26 
respectively for M = 240 and 320. The scales of the two ordinates are related by means 
of r* = 1.4 x 10~ 6 . We see from the figure that the field dependence of /zbd, particularly 
that of M = 240, well coincide with that of /i exp of DNA with the corresponding length we 
have specified, i.e., 4.3 and 6.5 kbp DNA for M = 240 and 320, respectively. This is quite 
satisfactory if we take into account the fact that only three adjustable parameters are used 
to draw the two sets of data in Fig. |3|. We therefore consider that the present BD model 
catches up even quantitatively essences of DNA dynamics under CFGE. 

When we compare /xbd with fi exp further in detail, however, there are some unsatisfactory 
features in the BD results. As compared with /z exp of 6.5 kbp DNA, an ^-independent branch 
of fiBD has not been ascertained in the weak field limit for the M = 320 chain. One possible 
reason for this may be the insufficient time of simulation (up to t = 10 6 in finite fields) for 
this length of chain in weak fields. Thus, within the limited CPU time, it is a rather hard 
task for the present BD simulation to judge the field dependence of /xbd; proportional to 
E or E 2 , in weak fields.@S0 In the strong field regime, on the other hand, /zbd tend to 



saturate, while, jj exp are still significantly increasing with E. These different tendencies are 
seen more clearly if we compare yU ex p measured in lOV/cm with the corresponding /xbd (not 
seen). The reason of this discrepancy is not clear at the moment. 

We have also examined /ibd of the chains in gels with a ge \ = 20 and 10. If only the volume 
ratio of gel is considered, they correspond roughly to 0.75% and 3% agaroses, respectively. 
Qualitatively, /zbd becomes larger and its independence becomes weaker for the larg er dge\ as 
expected, but these tendencies are quantitatively weaker than those observed experimentally. 
The origins of the discrepancy may be attributed to our model for gel, i.e., a perfectly rigid, 
regular jungle gym. 

IV. ELONGATION-CONTRACTION MOTION 

A. Motion of an individual chain 

In fields larger than a certain value depending on M and a ge i, chains in the present 3D 
BD simulation exhibit quasi-periodic behavior, i.e., they exhibit elongated and contracted 
shapes alternatively, as observed in the experiments0'lil'ii and in the previous 2D BD mod- 
els.SiH It is here emphasized that such quasi-periodic behavior has been observed within 
the field range where we have discussed the .E-dependence of /i in Fig. || We show a typical 
conformational change of a chain observed in one MD run of CFGE with the parameters 
M = 240, E = 0.032, and a ge i = 17 in Fig. g. Within a time window of the figure the 
elongation-contraction motion occurs three times at around t =2, 6 and 11 xlO 5 . At ranges 
between them centered around t — 4 and 8 xlO 5 a chain is in a rather compact form. In 
this range it can happen, though not so frequently, that the front and rear ends exchange 
their roles as seen around i = 9x 10 5 in the figure. It is also pointed out here that the front 
end of a chain moves in the field direction with almost a constant rate. 

B. Quasi-periodic behavior of Ri(t) and vc(t) 

In Fig. [5] we show the time evolution of the radius of longer principal axis, Ri(t), and the 
velocity of the center of mass, vo{t), in the MD run whose chain conformational change is 
shown in Fig. £2[ The fluctuation in Ri(t) is notable, showing successive A-shaped peaks. 
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We also note that just before the maxima of Ri(t), minima are observed in vc(t). These 
features are qualitatively in good agreement with those observed in the experiment. El 

The time evolution of Ri(t) and vo(t) in Fig. |5] looks quasi-periodic. Actually it was found 
by the experiment that the autocorrelation function of Ri(t), C^(t), exhibits a damped 
oscillation.!^ This is also the case for our present BD data. In Fig. |^ Crr(£) evaluated for 
various M are presented. For M > 80, Crr(£) exhibit the first undershoot below the line 
of Crr(^) = 0, and even the first overshoot is clearly seen for M > 160. From these results 
we may conclude that, in chains with M > 160 under E = 0.032, elongation-contraction 
motions occur quasi-periodically. In the inset of Fig. we show nearly linear dependence of 
the period of oscillation r = 4t on M, where t is defined as the time of the first intercept 
of Crr(£) with the abscissa!! 



C. Analyses on peak structures of Ri(t) 

Oana et argued that the steady-state time evolution of DNA under CFGE is clas- 
sified into two types of behavior. One is the elongation-contraction motion which looks 
apparently deterministic. We call the time branches in which such a motion is undergoing 
the deterministic ones (ranges centered at around t =2, 6 and 11 x 10 5 in Fig. The other 
is what is observed between two deterministic branches. We call them stochastic branches. 
Oana et al. approximately described the quasi-periodic behavior of Ri(t) in terms of a sim- 
plified function: peaks in Ri(t) are approximated by A-shape branches whose widths are 
assumed to obey a Gaussian distribution, and others by constants whose duration times are 
assumed to obey a Poisson distribution. They calculated Crr(^) by means of this model 
function for various cases and remarked that when XD becomes larger than unity a damped 
oscillation in Crr(£) becomes prominent. Here 1/A and D stand for the mean period of the 
stochastic branches and the mean width of the deterministic ones, respectively. 

Following the classification due to Oana et al, we analyze our BD results by making 
use of an algorithm which appropriately picks up peaks of Ri(t). For this purpose we first 
pick up time sets {t^m-i' *maxi ^min:r}' i- e -' the times of preceding local minimum of the i-th 
maximum, the i-th maximum itself, and the subsequent minimum. Then we select out peaks 
whose height relative to their adjacent minima is larger than a certain threshold value. The 
latter is chosen in such a way that the sum of periods of peak branches selected surmounts 
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more than 50% of the total duration of the observation. This procedure divides a whole time 
sequence of Ri(t) into the two types of branches whose numbers are equal to each other. It 
picks up ensemble of various peaks but not a limited numbers of peaks of special shapes. 
The number of peaks thus selected is more than a hundred for each set of parameters. 



D. Statistics of deterministic and stochastic branches 



For the set of times {t 7 ^^, t 7 ^^, t 7 ^^} of peaks picked up by the procedure described 
above we examine the distribution of peak widths (or periods of deterministic branch) D n = 
^min:r — ^min:i an d that of duration times between neighboring peaks (periods of stochastic 
branch) 1/A n = t^ in . T — £™7n r The results for M = 240, E = 0.032 and a ge \ = 17 are shown 
in Fig. 0. They numerically confirm the basic assumptions of Oana et ai, i.e., {D n } nearly 
obey a Gaussian distribution, while {1/A n } do a Poisson distribution. The mean and the 
variance of the former are given by D = 3.3 x 10 5 and A ~ 0.2D, respectively, while the 
mean width of the latter is given by XD ~ 1.7. The sum D + l/X should be equal to r = 4to 
evaluated from Crr^) in Sec. [IV Bj , which is verified within our numerical accuracy. The 
ratio XD ~ 1.6 fulfills the condition for a damped oscillation to be observed in Crr(£). 

The ratio XD is found to decrease below unity for M which is smaller than a certain 
value. In Fig. || D and 1/A are plotted against M for E = 0.032 and a gc \ = 17. The results 
strongly suggest that D oc M 1 and 1/A oc M° for large M. In Fig. |9|, the M-dependence of 
ratios XD and A/D are presented. As expected from the results shown in Fig. || XD is an 
approximately linearly increasing function for large M. On the other hand, A/D is almost 
constant (~ 0.3), indicating that the Gaussian distribution of {D n } is scaled by the mean 
D alone. It should be noted that, for M > 100 where XD > 1 holds, C^(t) clearly exhibits 
a damped oscillation as has been observed in Fig. |6|. Also the result 4t oc M for M > 100 
shown in the inset of Fig. ^| is in accordance with the fact that XD is larger than unity for 
such M. For sufficiently large M the deterministic branches dominate a whole sequence of 
the time evolution, and so 4t — D oc M holds. These results consistently indicate that there 
exists a crossover between the regimes with and without elongation-contraction motions at 
around M = 100 in CFGE with the parameters studied, i.e., E = 0.032 and a ge i = 17. 

We have also examined the field dependence of D and 1/A for a fixed M(= 240), the 



results of which are shown in Fig. 10. In strong fields E > 0.06, DE becomes to be saturated, 
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while E/X decreases with increasing E. As a consequence, DX is an increasing function of 
E, indicating that the deterministic branches, or in other words, the elongation-contraction 
motions survive and even dominate a whole evolution of the chain in sufficiently strong 
fields. 

Lastly, in Fig. [TT], the averaged time evolution of individual peaks of Ri(t) are plotted 
for various M and E. The abscissa and ordinate are normalized by the period r = 4to and 
the averaged peak height, respectively. It is seen that with these normalizations, the data 
with different M and E almost lie on a universal curve during the period from — 0.6r to 
0.4t. The results obtained above such as D oc M, E' 1 , A oc D(oc M) can be derived from 
the scaling behavior of Ri(t) in Fig. [H]. Such scaling behavior of Ri(t) including the fact 
that the normalized peak exhibits an anti-symmetric shape is just what was observed exper- 
imentally.El It was mentioned in Ref. 15 that the ratio of slope of the A-shape immediately 
after the maximum to that before it is approximately constant (~ 3) for all conditions they 
employed. The corresponding ratio of the BD result is ~ 2. The origin of this quantitative 
difference is not clear at the moment. 



V. DISCUSSIONS 



We have performed extensive simulations on the DNA constant-field gel electrophoresis 
(CFGE) by means of a new Brownian dynamics (BD) method which we have recently de- 
veloped. Our BD method is based on a generalized bead-spring model for a polymer in a 
3D space with rigid, regularly arrayed obstacles. By the word 'generalized' we mean that 
the 'spring force' is not linear but extremely nonlinear: it is infinite when the distance I 
between two beads connected by the 'spring' is larger than v2, while it is zero when / < y/2. 
Also in the model we strictly take into account of the excluded volume effects between any 
pair of beads whose diameter is put unity. These forces are supposed to be of a microscopic 
origin, i.e., molecular forces acting on monomers of a polymer. It should be emphasized that 
there involves no characteristic scale in these interactions except for the maximum distance 
between the connected beads (a/2) relative to the bead diameter (unity). 

In a vanishing field the viscosity £ (or At/( with At being the time unit of the BD 
dynamics) and the temperature T in eqn. [I] give rise to a characteristic scale of the model. 
The two parameters specify the viscous force and the random force from a solvent and 
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are related to each other by eqn. ||. As is discussed in the accompanying paper II, the 
characteristic scale is the 'persistence length' of the chain. Actually, with ( = 10, At = 1 
and A;bT B d = 2 the mean distance between two connected beads (cl* in eqn. ||) turns out to 
be about one six-th of the persistence length. Equating the latter to the persistence length 
of real DNA, the mean bead distance is roughly estimated as lOnm. Thus a bead in our BD 
model represents a potion of DNA consisting of a few tens of base-pairs. We may call our 
BD model with these values of parameters a bead-spring model of a semi-microscopic level. 
We will argue in paper III that the resultant chain dynamics, when it is coarse-grained in 
time and space, can be interpreted as dynamics of a charged, elastic string, and that the 
elasticity is due to the conformational entropy of the original, semi-microscopic chain. 

A value of the electric field, E = qE^ in eqn. |], is reasonably converted to that of the 
experimental CFGE on DNA as discussed in Sec. |TTT| . By the procedure adopted there the 
conversion of the field scale is essentially governed by that of the length scale as expressed 
in eqn. 10. The parameter c in eqn. |X0| has been adjusted for the chains with different 
sizes, but within the range that the results are compatible with the condition \pl > (I) > 1 
where (/) is the mean distance of connected beads. Then the ^-dependence of mobility /zbd 
has been shown to agree semi- quantitatively with that of /z exp once absolute magnitudes 
of the simulational and experimental /i's are adjusted at the vanishing field limit. This 
procedure fixes the conversion of the time scales of the simulation and the experiment. 
To our knowledge, such a semi- quantitative coincidence between the simulation and the 
experiment has not been so far achieved. 

We have also demonstrated that our BD simulation on CFGE in fields stronger than a 
certain crossover value reproduces the quasi-periodic evolution of a chain whose character- 
istic aspects are quite similar to those experimentally observed. Actually, as proposed by 
Oana et a/.El based on their experimental observations, a whole sequence of time evolution of 
the chain can be divided into two branches according to behavior of the longer radius of the 
chain Ri(t). One is the deterministic branch where Ri(t) exhibits a distinct peak structure 
which corresponds to the elongation-contraction motion of the chain, and the other is the 
stochastic branch between two neighboring deterministic branches where the chain exhibits 
a rather compact form. 

We have studied statistics of time durations of the two branches and have found that, as 
assumed by Oana et al, 111 those of the deterministic branch obey a Gaussian distribution 
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with mean D and variance A, while those of the stochastic branch do a Poisson distribution 
with the mean 1/A. Furthermore it is found that, in a fixed field, D oc M 1 , A oc D and 
1/A oc M° for large M. This means that the deterministic branches become more dominant 
in a whole time sequence, and that the quasi-periodicity in the chain dynamics becomes 
more distinct (see Fig. ^) for larger M. These results, combined with the observation that 
the chain moves to the field direction by the distance proportional to M in each deterministic 
branch, yield fi oc M°. This explains why CFGE cannot separate DNA according to their 
size L (oc M) for large L. We have also checked that, for a fixed M which is sufficiently 
large, the deterministic branches become dominant as E is increased at least within the field 
range of our BD simulation (see Fig. 10). 

By the present BD simulation a crossover between the regimes with and without 
elongation-contraction motions is expected to occur for a chain with M ~ 100 in E = 0.032 
and a ge i = 17. According to the length and field conversions described in Sec. |T|, this 
roughly corresponds to DNA of 2kbp in 1% agarose in E exp ~ 4V/cm. Experimently, on 
the other hand, the direct measurement of the elongation-contraction motion by the flu- 
orescent microscopy has been rather limited; the shortest is for T7 DNA (38kbp),S but 
mostly for T4 DNA (166kbp)0 and T2 DNA (164kbp)H But, there exist another experi- 
ment which strongly indicates the occurrence of the elongation-contraction motion, i.e., the 
antiresonance of mobility in the filed inversion gel electrophoresis. In such a experiment the 
antiresonance has been clearly observed for DNA with 9.42kbpS Therefore, our criterion 
mentioned above is considered to be a reasonable estimate for the elongation-contraction 
motion to occur. 

Lastly two further comments are in order. In the whole parameter ranges we have exam- 
ined, including the ones where the deterministic branches sufficiently dominate the stochas- 
tic ones, so-called 'hernias'-like configurations of a chain has been scarcely observed in the 
present BD simulation. This implies that they do not play a significant role in determining 
the saturation of fi at least for chains with moderate sizes we have studied. Another inter- 
esting observation has been already pointed out in Sec. [IV B| : the front end of a chain moves 
in the field direction with almost a constant rate (Fig. In relation to this, we note that 
the average of vc(t) over the deterministic branches and that over the stochastic branches 
are found to almost coincide with each other in a whole range of field we have examined. 
We shall interpret these results in paper III where the details of the chain dynamics in the 
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deterministic branches is discussed. 
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Figures 



FIG. 1: A snapshot of M = 10 chain simulated in the BD method. The spherical objects express 
beads with a radius 0.5. There is no overlap between any pair of beads. The bar-like objects are 
inserted only to visualize bonding of neighboring beads. 

FIG. 2: The scattering processes and the associated constraining forces to fulfill the constraint of 
eqn. |2| (a) and eqn. || (b). 

FIG. 3: Field dependence of mobilities. 

FIG. 4: A time evolution of projection of an M = 240 chain. E = 0.032 a ge \ = 17. Solid line 
stands for locus of % = 1, and dashed line for i = 240. 

FIG. 5: Time development of Ri(t) and vq under E = 0.032. 

FIG. 6: The autocorrelation functions Cnn(t). The dependence of period 4to on M is shown in 
the inset. 

FIG. 7: (a) Distributions of peak widths D n = in . r — in .j and (b) duration between neighboring 
peak minima 1/A„ = t^ :T - t" in:1 obtained for M = 240, E = 0.032, a gol = 17. 

FIG. 8: Dependence of D and 1/A on M. 

FIG. 9: Ratios XD and A/D. 

FIG. 10: Field dependence of DE, E/X and DX of M = 240 chain. 

FIG. 11: The averaged peak shape of Ri{t). The abscissa is normalized by the period r and the 
ordinate by the average of the peak height Ri(0). 
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FIG. 7 
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FIG. 9 
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FIG. 10 



20 
18 
16 
14 
12 
10 



M=240 






DE i" EH 






E/X -©-■■ 






DX : A : 




j ; 




































-o - 



10 



12 



11 

10 

9 



14 



£xl0 2 



29 



FIG. 11 
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